scenarios/scripts/followup/Prevalence plots_10.R

In statnet/PrEP-for-ASMM-and-adult-MSM: Network-Based Epidemic Modeling of HIV Transmission among MSM and ASMM populations

## whamp analysis file

library(EpiModelHIV)
library(EpiModelHPC)


s1 <- merge_simfiles(1, indir = "~/Campcl/scenarios/data/")
s2 <- merge_simfiles(2, indir = "~/Campcl/scenarios/data/")
s3 <- merge_simfiles(3, indir = "~/Campcl/scenarios/data/")
s4 <- merge_simfiles(4, indir = "~/Campcl/scenarios/data/")
s5 <- merge_simfiles(5, indir = "~/Campcl/scenarios/data/")
s6 <- merge_simfiles(6, indir = "~/Campcl/scenarios/data/")

# Figure 1 ----------------------------------------------------------------

# Line plots of cumulative PIA

steps <- 1:520

base.prev.all<-rep(NA,520)
base.prev.all_hi<-rep(NA,520)
base.prev.all_low<-rep(NA,520)

base.prev.msm<-rep(NA,520)
base.prev.msm_hi<-rep(NA,520)
base.prev.msm_low<-rep(NA,520)

base.prev.asmm<-rep(NA,520)
base.prev.asmm_hi<-rep(NA,520)
base.prev.asmm_low<-rep(NA,520)



##Intervention 1 adult only.
inter.1.prev.all<-rep(NA,520)
inter.1.prev.all_hi<-rep(NA,520)
inter.1.prev.all_low<-rep(NA,520)

inter.1.prev.msm<-rep(NA,520)
inter.1.prev.msm_hi<-rep(NA,520)
inter.1.prev.msm_low<-rep(NA,520)

inter.1.prev.asmm<-rep(NA,520)
inter.1.prev.asmm_hi<-rep(NA,520)
inter.1.prev.asmm_low<-rep(NA,520)


##Intervention 2 Adults plus ASMM @ 10%

inter.2.prev.all<-rep(NA,520)
inter.2.prev.all_hi<-rep(NA,520)
inter.2.prev.all_low<-rep(NA,520)

inter.2.prev.msm<-rep(NA,520)
inter.2.prev.msm_hi<-rep(NA,520)
inter.2.prev.msm_low<-rep(NA,520)

inter.2.prev.asmm<-rep(NA,520)
inter.2.prev.asmm_hi<-rep(NA,520)
inter.2.prev.asmm_low<-rep(NA,520)


##Intervention 3 Adults plus ASMM @ 20%

inter.3.prev.all<-rep(NA,520)
inter.3.prev.all_hi<-rep(NA,520)
inter.3.prev.all_low<-rep(NA,520)

inter.3.prev.msm<-rep(NA,520)
inter.3.prev.msm_hi<-rep(NA,520)
inter.3.prev.msm_low<-rep(NA,520)

inter.3.prev.asmm<-rep(NA,520)
inter.3.prev.asmm_hi<-rep(NA,520)
inter.3.prev.asmm_low<-rep(NA,520)


##Intervention 4 Adults plus ASMM @ 30%

inter.4.prev.all<-rep(NA,520)
inter.4.prev.all_hi<-rep(NA,520)
inter.4.prev.all_low<-rep(NA,520)

inter.4.prev.msm<-rep(NA,520)
inter.4.prev.msm_hi<-rep(NA,520)
inter.4.prev.msm_low<-rep(NA,520)

inter.4.prev.asmm<-rep(NA,520)
inter.4.prev.asmm_hi<-rep(NA,520)
inter.4.prev.asmm_low<-rep(NA,520)

##Intervention 5 Adults plus ASMM @ 40%

inter.5.prev.all<-rep(NA,520)
inter.5.prev.all_hi<-rep(NA,520)
inter.5.prev.all_low<-rep(NA,520)

inter.5.prev.msm<-rep(NA,520)
inter.5.prev.msm_hi<-rep(NA,520)
inter.5.prev.msm_low<-rep(NA,520)

inter.5.prev.asmm<-rep(NA,520)
inter.5.prev.asmm_hi<-rep(NA,520)
inter.5.prev.asmm_low<-rep(NA,520)


x.all <- x.1.all <- x.2.all <- x.3.all <- x.4.all <- x.5.all <- NULL
x.msm <- x.1.msm <- x.2.msm <- x.3.msm <- x.4.msm <- x.5.msm <- NULL 
x.asmm <- x.1.asmm <- x.2.asmm <- x.3.asmm <- x.4.asmm <- x.5.asmm  <- NULL

for (i in seq_along(steps)) {
  sim.base <- truncate_sim(s1, at = 8261)
  
  x.all<-sort(as.numeric(sim.base$epi$i.prev[i,1:100]))
  base.prev.all[i]<-mean(x.all)
  base.prev.all_hi[i]<-mean(x.all[97],x.all[98])
  base.prev.all_low[i]<-mean(x.all[3],x.all[2])
  
  x.msm<-sort(as.numeric(sim.base$epi$i.prev.msm[i,1:100]))
  base.prev.msm[i]<-mean(x.msm)
  base.prev.msm_hi[i]<-mean(x.msm[97],x.msm[98])
  base.prev.msm_low[i]<-mean(x.msm[3],x.msm[2])
  
  x.asmm<-sort(as.numeric(sim.base$epi$i.prev.age18[i,1:100]))
  base.prev.asmm[i]<-mean(x.asmm)
  base.prev.asmm_hi[i]<-mean(x.asmm[97],x.asmm[98])
  base.prev.asmm_low[i]<-mean(x.asmm[3],x.asmm[2])
  
##################.
  sim.int1 <- truncate_sim(s2, at = 8261)
  mn <- head(as.data.frame(sim.int1), steps[i])
  
  x.1.all<-sort(as.numeric(sim.int1$epi$i.prev[i,1:100]))
  inter.1.prev.all[i]<-mean(x.1.all)
  inter.1.prev.all_hi[i]<-mean(x.1.all[97],x.1.all[98])
  inter.1.prev.all_low[i]<-mean(x.1.all[3],x.1.all[2])
  
  x.1.msm<-sort(as.numeric(sim.int1$epi$i.prev.msm[i,1:100]))
  inter.1.prev.msm[i]<-mean(x.1.msm)
  inter.1.prev.msm_hi[i]<-mean(x.1.msm[97],x.1.msm[98])
  inter.1.prev.msm_low[i]<-mean(x.1.msm[3],x.1.msm[2])
  
  x.1.asmm<-sort(as.numeric(sim.int1$epi$i.prev.age18[i,1:100]))
  inter.1.prev.asmm[i]<-mean(x.1.asmm)
  inter.1.prev.asmm_hi[i]<-mean(x.1.asmm[97],x.1.asmm[98])
  inter.1.prev.asmm_low[i]<-mean(x.1.asmm[3],x.1.asmm[2])


  ##############.
  sim.int2 <- truncate_sim(s3, at = 8261)
  mn <- head(as.data.frame(sim.int2), steps[i])
  
  x.2.all<-sort(as.numeric(sim.int2$epi$i.prev[i,1:100]))
  inter.2.prev.all[i]<-mean(x.2.all)
  inter.2.prev.all_hi[i]<-mean(x.2.all[97],x.2.all[98])
  inter.2.prev.all_low[i]<-mean(x.2.all[3],x.2.all[2])
  
  x.2.msm<-sort(as.numeric(sim.int2$epi$i.prev.msm[i,1:100]))
  inter.2.prev.msm[i]<-mean(x.2.msm)
  inter.2.prev.msm_hi[i]<-mean(x.2.msm[97],x.2.msm[98])
  inter.2.prev.msm_low[i]<-mean(x.2.msm[3],x.2.msm[2])
  
  x.2.asmm<-sort(as.numeric(sim.int2$epi$i.prev.age18[i,1:100]))
  inter.2.prev.asmm[i]<-mean(x.2.asmm)
  inter.2.prev.asmm_hi[i]<-mean(x.2.asmm[97],x.2.asmm[98])
  inter.2.prev.asmm_low[i]<-mean(x.2.asmm[3],x.2.asmm[2])
  

  ##############.
  sim.int3 <- truncate_sim(s4, at = 8261)
  mn <- head(as.data.frame(sim.int3), steps[i])
  
  x.3.all<-sort(as.numeric(sim.int3$epi$i.prev[i,1:100]))
  inter.3.prev.all[i]<-mean(x.3.all)
  inter.3.prev.all_hi[i]<-mean(x.3.all[97],x.3.all[98])
  inter.3.prev.all_low[i]<-mean(x.3.all[3],x.3.all[2])
  
  x.3.msm<-sort(as.numeric(sim.int3$epi$i.prev.msm[i,1:100]))
  inter.3.prev.msm[i]<-mean(x.3.msm)
  inter.3.prev.msm_hi[i]<-mean(x.3.msm[97],x.3.msm[98])
  inter.3.prev.msm_low[i]<-mean(x.3.msm[3],x.3.msm[2])
  
  x.3.asmm<-sort(as.numeric(sim.int3$epi$i.prev.age18[i,1:100]))
  inter.3.prev.asmm[i]<-mean(x.3.asmm)
  inter.3.prev.asmm_hi[i]<-mean(x.3.asmm[97],x.3.asmm[98])
  inter.3.prev.asmm_low[i]<-mean(x.3.asmm[3],x.3.asmm[2])
  

  
##############.
  sim.int4 <- truncate_sim(s5, at = 8261)
  mn <- head(as.data.frame(sim.int4), steps[i])
  
  x.4.all<-sort(as.numeric(sim.int4$epi$i.prev[i,1:100]))
  inter.4.prev.all[i]<-mean(x.4.all)
  inter.4.prev.all_hi[i]<-mean(x.4.all[97],x.4.all[98])
  inter.4.prev.all_low[i]<-mean(x.4.all[3],x.4.all[2])
  
  x.4.msm<-sort(as.numeric(sim.int4$epi$i.prev.msm[i,1:100]))
  inter.4.prev.msm[i]<-mean(x.4.msm)
  inter.4.prev.msm_hi[i]<-mean(x.4.msm[97],x.4.msm[98])
  inter.4.prev.msm_low[i]<-mean(x.4.msm[3],x.4.msm[2])
  
  x.4.asmm<-sort(as.numeric(sim.int4$epi$i.prev.age18[i,1:100]))
  inter.4.prev.asmm[i]<-mean(x.4.asmm)
  inter.4.prev.asmm_hi[i]<-mean(x.4.asmm[97],x.4.asmm[98])
  inter.4.prev.asmm_low[i]<-mean(x.4.asmm[3],x.4.asmm[2])
  

}

base.prev.all <- base.prev.all * 100
base.prev.all_hi <- base.prev.all_hi * 100
base.prev.all_low <- base.prev.all_low * 100

base.prev.msm <- base.prev.msm * 100
base.prev.msm_hi <- base.prev.msm_hi * 100
base.prev.msm_low <- base.prev.msm_low * 100

base.prev.asmm <- base.prev.asmm * 100
base.prev.asmm_hi <- base.prev.asmm_hi * 100
base.prev.asmm_low <- base.prev.asmm_low * 100

############# 1. 
inter.1.prev.all <- inter.1.prev.all * 100
inter.1.prev.all_hi <- inter.1.prev.all_hi * 100
inter.1.prev.all_low <- inter.1.prev.all_low * 100

inter.1.prev.msm <- inter.1.prev.msm * 100
inter.1.prev.msm_hi <- inter.1.prev.msm_hi * 100
inter.1.prev.msm_low <- inter.1.prev.msm_low * 100

inter.1.prev.asmm <- inter.1.prev.asmm * 100
inter.1.prev.asmm_hi <- inter.1.prev.asmm_hi * 100
inter.1.prev.asmm_low <- inter.1.prev.asmm_low * 100

################# 2.
inter.2.prev.all <- inter.2.prev.all * 100
inter.2.prev.all_hi <- inter.2.prev.all_hi * 100
inter.2.prev.all_low <- inter.2.prev.all_low * 100

inter.2.prev.msm <- inter.2.prev.msm * 100
inter.2.prev.msm_hi <- inter.2.prev.msm_hi * 100
inter.2.prev.msm_low <- inter.2.prev.msm_low * 100

inter.2.prev.asmm <- inter.2.prev.asmm * 100
inter.2.prev.asmm_hi <- inter.2.prev.asmm_hi * 100
inter.2.prev.asmm_low <- inter.2.prev.asmm_low * 100

################# 3.
inter.3.prev.all <- inter.3.prev.all * 100
inter.3.prev.all_hi <- inter.3.prev.all_hi * 100
inter.3.prev.all_low <- inter.3.prev.all_low * 100

inter.3.prev.msm <- inter.3.prev.msm * 100
inter.3.prev.msm_hi <- inter.3.prev.msm_hi * 100
inter.3.prev.msm_low <- inter.3.prev.msm_low * 100

inter.3.prev.asmm <- inter.3.prev.asmm * 100
inter.3.prev.asmm_hi <- inter.3.prev.asmm_hi * 100
inter.3.prev.asmm_low <- inter.3.prev.asmm_low * 100

################# 4.
inter.4.prev.all <- inter.4.prev.all * 100
inter.4.prev.all_hi <- inter.4.prev.all_hi * 100
inter.4.prev.all_low <- inter.4.prev.all_low * 100

inter.4.prev.msm <- inter.4.prev.msm * 100
inter.4.prev.msm_hi <- inter.4.prev.msm_hi * 100
inter.4.prev.msm_low <- inter.4.prev.msm_low * 100

inter.4.prev.asmm <- inter.4.prev.asmm * 100
inter.4.prev.asmm_hi <- inter.4.prev.asmm_hi * 100
inter.4.prev.asmm_low <- inter.4.prev.asmm_low * 100


library(wesanderson)
pal <- wes_palette("Zissou")[5]
pal2<- wes_palette("Zissou")[2]
  
## For Paper
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_PrEP 40_0.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,40,5)
plot(base.prev.all, type = "n", ylim = c(0, 35), lwd = 3, col="black", axes=FALSE,
     main = "HIV Prevalence: Adult MSM PrEP - 40%, ASMM PrEP - 0%", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##ALL
xx <- c(1:(length(base.prev.all)), (length(base.prev.all)):1)
yy <- c(base.prev.all_low, rev(base.prev.all_hi))
polygon(xx, yy, col = EpiModel::transco("light blue", alpha = 0.5), border = NA)
lines(base.prev.all, lwd = 1, col = "light blue")

aa <- c(1:(length(inter.1.prev.all)), (length(inter.1.prev.all)):1)
bb <- c(inter.1.prev.all_low, rev(inter.1.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.3), border = NA)
lines(inter.1.prev.all, lwd = 2, col = "blue")

##MSM
xx <- c(1:(length(base.prev.msm)), (length(base.prev.msm)):1)
yy <- c(base.prev.msm_low, rev(base.prev.msm_hi))
polygon(xx, yy, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(base.prev.msm, lwd = 1, col = "red")

aa <- c(1:(length(inter.1.prev.msm)), (length(inter.1.prev.msm)):1)
bb <- c(inter.1.prev.msm_low, rev(inter.1.prev.msm_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(inter.1.prev.msm, lwd = 1, col = "red")

##ASMM
xx <- c(1:(length(base.prev.asmm)), (length(base.prev.asmm)):1)
yy <- c(base.prev.asmm_low, rev(base.prev.asmm_hi))
polygon(xx, yy, col = EpiModel::transco("light green", alpha = 0.3), border = NA)
lines(base.prev.asmm, lwd = 1, col = "light green")

aa <- c(1:(length(inter.1.prev.asmm)), (length(inter.1.prev.asmm)):1)
bb <- c(inter.1.prev.asmm_low, rev(inter.1.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.3), border = NA)
lines(inter.1.prev.asmm, lwd = 2, col = "green")

legend(x=360,y=35, c("Total Population", "Adult MSM", "ASMM age 18"),text.col=c("blue","red", "green"), 
       border = "black", lty=1, col=c("blue","red", "green"))


dev.off()

###############################################################################################
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_PrEP 40_10.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,40,5)
plot(base.prev.all, type = "n", ylim = c(0, 35), lwd = 3, col="black", axes=FALSE,
     main = "HIV Prevalence: Adult MSM PrEP - 40%, ASMM PrEP - 10%", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##ALL
xx <- c(1:(length(base.prev.all)), (length(base.prev.all)):1)
yy <- c(base.prev.all_low, rev(base.prev.all_hi))
polygon(xx, yy, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(base.prev.all, lwd = 1, col = "red")

aa <- c(1:(length(inter.2.prev.all)), (length(inter.2.prev.all)):1)
bb <- c(inter.2.prev.all_low, rev(inter.2.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(inter.2.prev.all, lwd = 2, col = "red")

##MSM
xx <- c(1:(length(base.prev.msm)), (length(base.prev.msm)):1)
yy <- c(base.prev.msm_low, rev(base.prev.msm_hi))
polygon(xx, yy, col = EpiModel::transco("light blue", alpha = 0.3), border = NA)
lines(base.prev.msm, lwd = 2, col = "light blue")

aa <- c(1:(length(inter.2.prev.msm)), (length(inter.2.prev.msm)):1)
bb <- c(inter.2.prev.msm_low, rev(inter.2.prev.msm_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.3), border = NA)
lines(inter.2.prev.msm, lwd = 2, col = "blue")

##ASMM
xx <- c(1:(length(base.prev.asmm)), (length(base.prev.asmm)):1)
yy <- c(base.prev.asmm_low, rev(base.prev.asmm_hi))
polygon(xx, yy, col = EpiModel::transco("light green", alpha = 0.3), border = NA)
lines(base.prev.asmm, lwd = 2, col = "light green")

aa <- c(1:(length(inter.2.prev.asmm)), (length(inter.2.prev.asmm)):1)
bb <- c(inter.2.prev.asmm_low, rev(inter.2.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.3), border = NA)
lines(inter.2.prev.asmm, lwd = 2, col = "green")

legend(x=360,y=35, c("Total Population", "Adult MSM", "ASMM age 18"),text.col=c("blue","red", "green"), 
       border = "black", lty=1, col=c("blue","red", "green"))



dev.off()

###############################################################################################
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_PrEP 40_20.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,40,5)
plot(base.prev.all, type = "n", ylim = c(0, 35), lwd = 3, col="black", axes=FALSE,
     main = "HIV Prevalence: Adult MSM PrEP - 40%, ASMM PrEP - 20%", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##ALL
xx <- c(1:(length(base.prev.all)), (length(base.prev.all)):1)
yy <- c(base.prev.all_low, rev(base.prev.all_hi))
polygon(xx, yy, col = EpiModel::transco("light blue", alpha = 0.3), border = NA)
lines(base.prev.all, lwd = 2, col = "light blue")

aa <- c(1:(length(inter.3.prev.all)), (length(inter.3.prev.all)):1)
bb <- c(inter.3.prev.all_low, rev(inter.3.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.3), border = NA)
lines(inter.3.prev.all, lwd = 2, col = "blue")

##MSM
xx <- c(1:(length(base.prev.msm)), (length(base.prev.msm)):1)
yy <- c(base.prev.msm_low, rev(base.prev.msm_hi))
polygon(xx, yy, col = EpiModel::transco("pink", alpha = 0.3), border = NA)
lines(base.prev.msm, lwd = 2, col = "pink")

aa <- c(1:(length(inter.3.prev.msm)), (length(inter.3.prev.msm)):1)
bb <- c(inter.3.prev.msm_low, rev(inter.3.prev.msm_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(inter.3.prev.msm, lwd = 2, col = "red")

##ASMM
xx <- c(1:(length(base.prev.asmm)), (length(base.prev.asmm)):1)
yy <- c(base.prev.asmm_low, rev(base.prev.asmm_hi))
polygon(xx, yy, col = EpiModel::transco("light green", alpha = 0.3), border = NA)
lines(base.prev.asmm, lwd = 2, col = "light green")

aa <- c(1:(length(inter.3.prev.asmm)), (length(inter.3.prev.asmm)):1)
bb <- c(inter.3.prev.asmm_low, rev(inter.3.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.3), border = NA)
lines(inter.3.prev.asmm, lwd = 2, col = "green")

legend(x=360,y=35, c("Total Population", "Adult MSM", "ASMM age 18"),text.col=c("blue","red", "green"), 
       border = "black", lty=1, col=c("blue","red", "green"))


dev.off()

###############################################################################################
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_PrEP 40_30.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,40,5)
plot(base.prev.all, type = "n", ylim = c(0, 35), lwd = 3, col="black", axes=FALSE,
     main = "HIV Prevalence: Adult MSM PrEP - 40%, ASMM PrEP - 30%", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##ALL
xx <- c(1:(length(base.prev.all)), (length(base.prev.all)):1)
yy <- c(base.prev.all_low, rev(base.prev.all_hi))
polygon(xx, yy, col = EpiModel::transco("light blue", alpha = 0.3), border = NA)
lines(base.prev.all, lwd = 2, col = "light blue")

aa <- c(1:(length(inter.4.prev.all)), (length(inter.4.prev.all)):1)
bb <- c(inter.4.prev.all_low, rev(inter.4.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.3), border = NA)
lines(inter.4.prev.all, lwd = 2, col = "blue")

##MSM
xx <- c(1:(length(base.prev.msm)), (length(base.prev.msm)):1)
yy <- c(base.prev.msm_low, rev(base.prev.msm_hi))
polygon(xx, yy, col = EpiModel::transco("pink", alpha = 0.3), border = NA)
lines(base.prev.msm, lwd = 2, col = "pink")

aa <- c(1:(length(inter.4.prev.msm)), (length(inter.4.prev.msm)):1)
bb <- c(inter.4.prev.msm_low, rev(inter.4.prev.msm_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(inter.4.prev.msm, lwd = 2, col = "red")

##ASMM
xx <- c(1:(length(base.prev.asmm)), (length(base.prev.asmm)):1)
yy <- c(base.prev.asmm_low, rev(base.prev.asmm_hi))
polygon(xx, yy, col = EpiModel::transco("light green", alpha = 0.3), border = NA)
lines(base.prev.asmm, lwd = 2, col = "light green")

aa <- c(1:(length(inter.4.prev.asmm)), (length(inter.4.prev.asmm)):1)
bb <- c(inter.4.prev.asmm_low, rev(inter.4.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.3), border = NA)
lines(inter.4.prev.asmm, lwd = 2, col = "green")

legend(x=360,y=35, c("Total Population", "Adult MSM", "ASMM age 18"),text.col=c("blue","red", "green"), 
       border = "black", lty=1, col=c("blue","red", "green"))


dev.off()


###############################################################################################
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_Prev all.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,40,5)
plot(base.prev.all, type = "n", ylim = c(0, 35), lwd = 3, col="black", axes=FALSE,
     main = "Population HIV Prevalence", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##Baseline
xx <- c(1:(length(base.prev.all)), (length(base.prev.all)):1)
yy <- c(base.prev.all_low, rev(base.prev.all_hi))
polygon(xx, yy, col = EpiModel::transco("grey", alpha = 0.3), border = NA)
lines(base.prev.all, lwd = 2, col = "black")

#INT1.
aa <- c(1:(length(inter.1.prev.all)), (length(inter.1.prev.all)):1)
bb <- c(inter.1.prev.all_low, rev(inter.1.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(inter.1.prev.all, lwd = 2, col = "red")

#INT2.
aa <- c(1:(length(inter.2.prev.all)), (length(inter.2.prev.all)):1)
bb <- c(inter.2.prev.all_low, rev(inter.2.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("orange", alpha = 0.3), border = NA)
lines(inter.2.prev.all, lwd = 2, col = "orange")

#INT3.
aa <- c(1:(length(inter.3.prev.all)), (length(inter.3.prev.all)):1)
bb <- c(inter.3.prev.all_low, rev(inter.3.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.3), border = NA)
lines(inter.3.prev.all, lwd = 2, col = "green")

#INT4.
aa <- c(1:(length(inter.4.prev.all)), (length(inter.4.prev.all)):1)
bb <- c(inter.4.prev.all_low, rev(inter.4.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.3), border = NA)
lines(inter.4.prev.all, lwd = 2, col = "blue")

#INT5.
#aa <- c(1:(length(inter.5.prev.all)), (length(inter.5.prev.all)):1)
#bb <- c(inter.5.prev.all_low, rev(inter.5.prev.all_hi))
#polygon(aa, bb, col = EpiModel::transco("violet", alpha = 0.3), border = NA)
#lines(inter.5.prev.all, lwd = 2, col = "violet")

legend(x=360,y=35, c("No PrEP", "Adult MSM PrEP 40%", "+ ASMM PrEP 10%", "+ ASMM PrEP 20%", "+ ASMM PrEP 30%"),text.col=c("black","red", "orange", "green", "blue"), 
       border = "black", lty=1, col=c("black","red", "orange", "green", "blue"))

dev.off()



###############################################################################################
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_Prev all 2 int.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,40,5)
plot(base.prev.all, type = "n", ylim = c(0, 35), lwd = 3, col="black", axes=FALSE,
     main = "HIV Prevalence with no PrEP, adult MSM PrEP and ASMM PrEP", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##Baseline
xx <- c(1:(length(base.prev.all)), (length(base.prev.all)):1)
yy <- c(base.prev.all_low, rev(base.prev.all_hi))
polygon(xx, yy, col = EpiModel::transco("grey", alpha = 0.3), border = NA)
lines(base.prev.all, lwd = 2, col = "black")

#INT1.
aa <- c(1:(length(inter.1.prev.all)), (length(inter.1.prev.all)):1)
bb <- c(inter.1.prev.all_low, rev(inter.1.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.3), border = NA)
lines(inter.1.prev.all, lwd = 2, col = "red")

#INT2.
#aa <- c(1:(length(inter.2.prev.all)), (length(inter.2.prev.all)):1)
#bb <- c(inter.2.prev.all_low, rev(inter.2.prev.all_hi))
#polygon(aa, bb, col = EpiModel::transco("orange", alpha = 0.3), border = NA)
#lines(inter.2.prev.all, lwd = 2, col = "orange")

#INT3.
#aa <- c(1:(length(inter.3.prev.all)), (length(inter.3.prev.all)):1)
#bb <- c(inter.3.prev.all_low, rev(inter.3.prev.all_hi))
#polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.3), border = NA)
#lines(inter.3.prev.all, lwd = 2, col = "green")

#INT4.
aa <- c(1:(length(inter.4.prev.all)), (length(inter.4.prev.all)):1)
bb <- c(inter.4.prev.all_low, rev(inter.4.prev.all_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.3), border = NA)
lines(inter.4.prev.all, lwd = 2, col = "blue")

#INT5.
#aa <- c(1:(length(inter.5.prev.all)), (length(inter.5.prev.all)):1)
#bb <- c(inter.5.prev.all_low, rev(inter.5.prev.all_hi))
#polygon(aa, bb, col = EpiModel::transco("violet", alpha = 0.3), border = NA)
#lines(inter.5.prev.all, lwd = 2, col = "violet")

legend(x=360,y=35, c("No PrEP", "Adult MSM PrEP 40%", "+ ASMM PrEP 30%"),text.col=c("black","red", "blue"), 
       border = "black", lty=1, col=c("black","red", "blue"))


dev.off()


###############################################################################################
#pdf(file = "Fig1.pdf", height = 6, width = 12, pointsize = 16)
tiff(filename = "Fig_Prev ASMM.tiff", height = 7, width = 9, units = "in", res = 250)

par(mfrow = c(1,1), mar = c(3,3,2.5,5), mgp = c(2,1,0))
xticks <- seq(0, 520, 52)
yticks <- seq(0,10,1)
plot(base.prev.asmm, type = "n", ylim = c(0, 10), lwd = 3, col="black", axes=FALSE,
     main = "HIV prevalence among ASMM age 18 under five PrEP interventions", xlab = "Years", ylab = "HIV Prevalence")
axis(1, at = xticks, labels = c("0","1","2","3","4","5","6","7","8","9","10"), col.axis="black", las=1, cex.axis=.6)
axis(2, at = yticks, col.axis="black", las=1, cex.axis=.6)

##Baseline
xx <- c(1:(length(base.prev.asmm)), (length(base.prev.asmm)):1)
yy <- c(base.prev.asmm_low, rev(base.prev.asmm_hi))
polygon(xx, yy, col = EpiModel::transco("grey", alpha = 0.2), border = NA)
lines(base.prev.asmm, lwd = 2, col = "black")

#INT1.
aa <- c(1:(length(inter.1.prev.asmm)), (length(inter.1.prev.asmm)):1)
bb <- c(inter.1.prev.asmm_low, rev(inter.1.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("red", alpha = 0.2), border = NA)
lines(inter.1.prev.asmm, lwd = 2, col = "red")

#INT2.
aa <- c(1:(length(inter.2.prev.asmm)), (length(inter.2.prev.asmm)):1)
bb <- c(inter.2.prev.asmm_low, rev(inter.2.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("orange", alpha = 0.2), border = NA)
lines(inter.2.prev.asmm, lwd = 2, col = "orange")

#INT3.
aa <- c(1:(length(inter.3.prev.asmm)), (length(inter.3.prev.asmm)):1)
bb <- c(inter.3.prev.asmm_low, rev(inter.3.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("green", alpha = 0.2), border = NA)
lines(inter.3.prev.asmm, lwd = 2, col = "green")

#INT4.
aa <- c(1:(length(inter.4.prev.asmm)), (length(inter.4.prev.asmm)):1)
bb <- c(inter.4.prev.asmm_low, rev(inter.4.prev.asmm_hi))
polygon(aa, bb, col = EpiModel::transco("blue", alpha = 0.2), border = NA)
lines(inter.4.prev.asmm, lwd = 2, col = "blue")

#INT5.
#aa <- c(1:(length(inter.5.prev.all)), (length(inter.5.prev.all)):1)
#bb <- c(inter.5.prev.all_low, rev(inter.5.prev.all_hi))
#polygon(aa, bb, col = EpiModel::transco("violet", alpha = 0.3), border = NA)
#lines(inter.5.prev.all, lwd = 2, col = "violet")

legend(x=360,y=10, c("No PrEP", "Adult MSM PrEP 40%", "+ ASMM PrEP 10%", "+ ASMM PrEP 20%", "+ ASMM PrEP 30%"),text.col=c("black","red", "orange", "green", "blue"), 
       border = "black", lty=1, col=c("black","red", "orange", "green", "blue"))


dev.off()

##Prevalence Table

base.prev.all <- tail(base.prev.all,1)
base.prev.all_low <- tail(base.prev.all_low,1)
base.prev.all_hi <- tail(base.prev.all_hi,1)

############# 1. 
inter.1.prev.all <- tail(inter.1.prev.all,1)
inter.1.prev.all_low <- tail(inter.1.prev.all_low,1)
inter.1.prev.all_hi <- tail(inter.1.prev.all_hi,1)

############# 2. 
inter.2.prev.all <- tail(inter.2.prev.all,1)
inter.2.prev.all_low <- tail(inter.2.prev.all_low,1)
inter.2.prev.all_hi <- tail(inter.2.prev.all_hi,1)

############# 3. 
inter.3.prev.all <- tail(inter.3.prev.all,1)
inter.3.prev.all_low <- tail(inter.3.prev.all_low,1)
inter.3.prev.all_hi <- tail(inter.3.prev.all_hi,1)


############# 4. 
inter.4.prev.all <- tail(inter.4.prev.all,1)
inter.4.prev.all_low <- tail(inter.4.prev.all_low,1)
inter.4.prev.all_hi <- tail(inter.4.prev.all_hi,1)



Prevalence.table.10yr <-cbind(
base.prev.all,
base.prev.all_low,
base.prev.all_hi,
inter.1.prev.all,
inter.1.prev.all_low,
inter.1.prev.all_hi,
inter.2.prev.all,
inter.2.prev.all_low,
inter.2.prev.all_hi,
inter.3.prev.all,
inter.3.prev.all_low,
inter.3.prev.all_hi,
inter.4.prev.all,
inter.4.prev.all_low,
inter.4.prev.all_hi)



library(xlsx) #load the package

write.xlsx(x = Prevalence.table.10yr, file = "/homes/dth2/Campcl/scenarios/out/Prevalence_table_10yr.xlsx",
           sheetName = "output", row.names = FALSE)



##Prevalence table for ASMM.


base.prev.asmm <- tail(base.prev.asmm,1)
base.prev.asmm_low <- tail(base.prev.asmm_low,1)
base.prev.asmm_hi <- tail(base.prev.asmm_hi,1)

############# 1. 
inter.1.prev.asmm <- tail(inter.1.prev.asmm,1)
inter.1.prev.asmm_low <- tail(inter.1.prev.asmm_low,1)
inter.1.prev.asmm_hi <- tail(inter.1.prev.asmm_hi,1)

############# 2. 
inter.2.prev.asmm <- tail(inter.2.prev.asmm,1)
inter.2.prev.asmm_low <- tail(inter.2.prev.asmm_low,1)
inter.2.prev.asmm_hi <- tail(inter.2.prev.asmm_hi,1)

############# 3. 
inter.3.prev.asmm <- tail(inter.3.prev.asmm,1)
inter.3.prev.asmm_low <- tail(inter.3.prev.asmm_low,1)
inter.3.prev.asmm_hi <- tail(inter.3.prev.asmm_hi,1)


############# 4. 
inter.4.prev.asmm <- tail(inter.4.prev.asmm,1)
inter.4.prev.asmm_low <- tail(inter.4.prev.asmm_low,1)
inter.4.prev.asmm_hi <- tail(inter.4.prev.asmm_hi,1)



Prevalence.asmm.table.10yr <-cbind(
  base.prev.asmm,
  base.prev.asmm_low,
  base.prev.asmm_hi,
  inter.1.prev.asmm,
  inter.1.prev.asmm_low,
  inter.1.prev.asmm_hi,
  inter.2.prev.asmm,
  inter.2.prev.asmm_low,
  inter.2.prev.asmm_hi,
  inter.3.prev.asmm,
  inter.3.prev.asmm_low,
  inter.3.prev.asmm_hi,
  inter.4.prev.asmm,
  inter.4.prev.asmm_low,
  inter.4.prev.asmm_hi)



library(xlsx) #load the package

write.xlsx(x = Prevalence.asmm.table.10yr, file = "/homes/dth2/Campcl/scenarios/out/Prevalence_asmm_table_ASMM_10yr.xlsx",
           sheetName = "output", row.names = FALSE)

Prevalence.asmm.table.10yr